Optimal series representations for numerical path integral simulations 
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By means of the Ito-Nisio theorem, we introduce and discuss a general approach to series repre- 
sentations of path integrals. We then argue that the optimal basis for both "primitive" and partial 
averaged approaches is the Wiener sine-Fourier basis. The present analysis also suggests a new 
approach to improving the convergence of primitive path integral methods. Current work indi- 
cates that this new technique, the "reweighted" method, converges as the cube of the number of 
path variables for "smooth" potentials. The technique is based on a special way of approximating 
the Brownian bridge which enters the Feynman-Kag formula and it does not require the Gaussian 
transform of the potential for its implementation. 
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I. INTRODUCTION 

Numerical simulations based on the path integral ap- 
proach have proved highly successful in the calculation 
of thermodynamic properties for complex, many-body 
quantum systems (see Refs. 0j| and tke cited bibjliog- 
raphy). Mainly the result of FeynmarJj and Kagjfl the 
centerpiece of the theory is the fact that the density 
matrix of a monodimensional system can be written as 
the expectation value of a suitable functional of a stan- 
dard Brownian bridge {B^, < u < 1}. More precisely, 
if {Bu, u > 0} is a standard Brownian motion start- 
ing at zero, then the Brownian bridge is the stochastic 
process {Bu\ Bi = 0, < w <. 1} i.e., a Brownian mo- 
tion conditioned on Bi — O.u In this paper, we shall 
reserve the symbol E to denote the expected value (av- 
erage value) of a certain random variable against the un- 
derlying probability measure of the Brownian bridge B^. 
For a monodimensional canonical ensemble character- 
ized by the inverse temperature (3 = l/(fcsT) and made 
up of identical particles of mass mo moving in the po- 
tentia^^a:), the Feynman-Kag density matrix formula 



Pfp{x,x';l3) 



= E exp ■ 



lo 



xo{u) + 




(1) 

where pfp{x,x';(3) stands for the density matrix of a 
similar free particle canonical ensemble, while Xf){u) is 
a shorthand for x + {x' — x)u. 

Current research is focused on the development of ac- 
curate, finite-dimensional approximations of the stochas- 
tic integrals that appear in Eq. 1 and in related thermo- 
dynamic expressions. The importance of Eq. 1 as given 
here consists of the fact that the Brownian motion, hence 
the Brownian bridge, are well understood mathematical 
objects, which can be simulated by a variety of means. 
The discussion in the present paper is based on the ran- 
dom series technique as a general representation scheme 
for the Brownian bridge B^. The approach is particularly 
interesting, as it is directly related to the "path integral" 



concept, |-and can be justified by means of the Ito-Nisio 
theoremJH whose statement is presented in Appendix A. 

We consider a number of questions related to the ran- 
dom series implementation of the Feynman-Kag formula. 
The so called primitives and partial averagingpl tech- 
niques, developed initially for the Fourier path integral 
(FPI) method,El are generalized here for arbitrary se- 
ries representations. Then, we address the question of 
whether or not there exists a preferred basis within which 
to implement the two techniques. We present strong ev- 
idence suggesting that the fastest convergent series for 
each method is the Wiener series on which the Fourier 
path integral approach is based. Finally, we introduce a 
new, non-averaging technique called the reweighted FPI 
method in order to improve the convergence of primitive 
FPI. 

Motivated by the optimality of the Wiener series, 
we undertake the task of establishing numerically the 
asymptotic rate of convergence for the three FPI meth- 
ods: the primitive FPI, the partial averaging FPI (PA- 
FPI), and the reweighted FPI (RW-FPI). The asymptotic 
rate of lUpavergence of the primitive FPI was extensively 
studiedclliy and is known to be 0{l/n) for sufficiently 
smooth potentials. However, there are at present no an- 
alytical or numerical studies concerning the exact asymp- 
totic behavior of the PA-FPI method. For the particu- 
lar case of the harmonic oscillator, it is known that the 
asymptotic rate of convergence is 0{l/ii?). (The reader 
should not mistake the full PA-FPI for the so called gra- 
dient corrected PA-FPI, which was shown to converge as 
fast as ©(l/n^) in Ref. |l^ for potentials having continu- 
ous second-order derivatives). To cope with the numer- 
ical difficulties encountered, we develop a Monte Carlo 
technique which allows us to study the asymptotic be- 
havior of the PA-FPI and RW-FPI methods, at least for 
single-well potentials. With its help, we find strong nu- 
merical evidence suggesting that the asymptotic rate of 
convergence for both PA-FPI and RW-FPI approaches is 
O^i/rc') for sufficiently smooth potentials. To our knowl- 
edge, RW-FPI thus becomes the most rapidly convergent 
method among those that leave the original potential un- 
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changed. 

The error analysis performed in Appendix E allows us 
to introduce what we call "accelerated" estimators, which 
are capable of improving the rate of convergence of any of 
the aforementioned methods from 0{l/n°') to 0{l/n"^^) 
for the first-order correction, and to 0{l/n"^'^) for the 
second-order correction, respectively. Although there is a 
price paid in the form of an increase in the variance of the 
respective estimators, the first-order correction appears 
suitable for general applications. 



II. SERIES REPRESENTATIONS OF THE 
BROWNIAN BRIDGE 

The most general series representation of the Brown- 
ian bridge is given by the Ito-Nisio theorem, the explicit 
statement of which is presented in Appendix A. We begin 
by assuming that we are given {Afe(T)}j.>i, a system of 
functions on the interval [0, 1], which, together with the 
constant function, Ao(t) = 1, make up an orthonormal 
basis in L^[0, 1]. If is the space of infinite sequences 
a = (ai, 02, . . .) and 



P[d] = II fi{ak) 



(2) 



fe=i 



is the (unique) probability measure on 17 such that the 
coordinate maps a are independent identically dis- 

tributed variables with distribution probability 



H{ak e A) = 



1 



-'/2dz 



then, 



S2(a) = VafcAfc(u), 0<ii<l 



k=l 



(3) 



(4) 



i.e., the right-hand side random series is equal in distri- 
bution to a standard Brownian bridge. Therefore, the 
notation B^{d) in (Q) is appropriate and allows us to in- 
terpret the Brownian bridge as a collection of random 
functions of argument a, indexed by u. 

Using the Ito-Nisio representation of the Brownian 
bridge, the Feynman-Kag formula (|^) takes the form 



pix,x';/3) 
Pfp{x,x';/3) 



dP[a]exp<^ - 13 V xo{u) + 



> akAk(u) 

t-[ 



(5) 



To reinforce the formula (^, consider the func- 
tions {'\/2cos(fc7rr)}fe>i, which, together with the con- 
stant function, make up a complete orthonormal system 
ofL2[o, 1]. Since 



-\/2cos(A:7rr)dr = 



2 sin(fc7rM) 



the Ito-Nisio theorem implies that 

sin(fc7ru) 



B 



k=l 



< u < 1, 



(6) 



so that the Feynman-Kag formula becomes 
p{x,x';f]) 



Pfp{x,x'](3) 



dP\ 



a\ exp 



(3 I V 

10 



xo(u) + 



afcCTfc sin(fc7ni) 



fe=l 



Au 



(7) 



where 



2 _ 2[ih^ 1 

— 9 "TT ■ 

Equation (7), derived here as a special case of the 
Ito-Nisio theorem, is the so-called Fourier path inte- 
gral method.Ei Historically, the sine-Fourier representa- 
tion was one of tiie first explicit constructions of the 
Brownian motion.O Following the mathematical litera- 
ture, we shall call it the Wiener construction after the 
name of its author, even though the original FPI method 
was deduced using arguments other than those presented 
here. 

The "primitive" series representation method consists 
of approximating the Brownian bridge by the n-th order 
partial sum of the series (0). Thus, 



Pfp{x,x';(3) 



dP[a] exp 



/3 / V 
Jo 



xoiu) 



Too 



fc=i 



afeAfe(M) 



(8) 



An immediate question arises: What is the best choice 
of functions Ai(u), i > 1, independent of potential, 
such that (H) has the fastest convergence? Although the 
phrase "independent of potential" carries ambiguities, in 
the remainder of this section we shall provide a more 
precise statement of the problem. 

We start with the observation that the Wiener basis is 
the only basis for which both Xi{u), and their primitives 
Ai(u), i > 1 are orthogonal. Indeed, let us notice that 
by construction, Ai(0) = for « > 1 and that 



A.(l) 



Ai(r)Ao(r)dT = Vi>l 



by orthogonality and the fact that Ao (t) = 1 . The unique 
basis, Ai(u), for which 



A,(T)A,(T)dr-0, Vz^j 



Ai(r)Aj(r)dr = V i, j > 1 

A,(0) = A,(l) =0, Vz>l 
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is made (up to a multiplication factor) of the eigenfunc- 
tions of the Dirichlct problem: 



1 



AA,{u) = e,A,(M), A,(0) = A,(l) = 0, 



as follows from the associated Dirichlet variational prin- 
ciple and the non-degeneracy of the spectrum of the "par- 
ticle in a box problem" . But that basis is precisely the 
Wiener basis. 

The orthogonality of the primitives, Xi{u), suggests 
that the Wiener basis is (in a sense that will be made 
clear below) optimal for the representation of the Brow- 
nian bridge. Let us define 



^::(a) = ^afeAfe(w) and K(a) 



k=l 



E 

k=n+l 



afcAfe(u) 



as the n-th order partial sum in (^ and the correspond- 
ing "tail" series, respectively. In terms of these sums, the 
Brownian bridge is expressed as B^{a) = 5" (a) -|- B!^{a). 
Obviously, and are independent. Moreover, a 
standard theorem regarding the sum of independent 
Gaussian distributed random variables shows that 
and S"" are again Gaussian distributed random variables 
of mean zero and variances 

oo n 
k=n+l k=l 

respectively. By independence, we have the equality 

E{By = e{b:)^ + E{s:r - ^(i - u), (9) 



where we used the fact that the variance of the Brownian 
bridge does not depend upon the series representation 
and so, it can be computed by using any convenient basis 
(e.g. the sine-Fourier basis). 

A natural way of measuring the quality of the approx- 
imation ^"(a) w B^l{a) is the value of the time average 
of the variances of the tails 



Jo 



s:ydu^ / EiB:rdu = 

Jo 

n 

u{l — u) — Afc(?i)^ du. 

k=l 



(10) 



Intuitively, the best approximating series is the one that 
minimizes the functional (10) for each n (we shall show 
that the answer is indeed a series) . More clearly, we want 
to find {Afe(r); k € l,n}, the system of functions on the 
interval [0, 1] which, together with the constant function 
Ao(t) — 1, make up an orthonormal system in L^[0, 1] 
and which realizes the maximum of the functional 



G(Ai,, 



n „i 

.^«) = E / 

k=l -^0 



AUu)du. 



(11) 



Since the system {\/2 cos(fc7rT)}fc>i together with the 
constant function make up a complete orthonormal sys- 
tem of L^[0, 1], we may write 



^ \/2cos(/7rM) / AA;(T)-\/2cos(/7rr)dr. 
,1 Jo 



Afe(u) 



Replacing this in (|1 1|) , we obtain 



G(Ai 



A„) = E/ A^Hd--EEE2/ / 

fc=l-^0 k=ll=l 7 = 1 -^0 



Xk{T)\ki0) cos(Z7rr) cos(j7r0)dTd0 x 



^ 2 sin(^7rw) sin(j7ru) 



k=l 



1=1 



2 cos(/7rT) cos(/7r6') 



drd6', 



(12) 



where we used the fact that the system {-\/2 sin(fc7rr)}fc>i 
is also orthonormal. From the theory of integral equa- 
tions with symmetric kernels, we learn that the maxi- 
mum of (12) is realized on the set of the n eigenfunctions 
having the largest eigenvalues. Since the kernel is al- 
ready in the series representation form, the maximum 
of our problem is J2k=i ^/i'^k)^ and is attained on the 
(orthonormal) functions 

Afc('u) = cos(fc7rM) fc G l,n. 

It follows that the Wiener representation is the unique 
series for which the time-average of the variance of the 



tail series reaches the minimum value of 



.1 n 
Jo 6 ^ TT^) 



fc2' 



(13) 



However, there is a direct connection between the 
asymptotic rate of convergence of the primitive method 
and the quantity E(i?" — S"")^, a connection that is given 
by formula (|2^) and is analyzed in Section III A. It allows 
us to conclude that the Wiener representation is the best 
series for general use in the primitive method. 
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III. IMPROVEMENTS IN THE PRIMITIVE 
FOURIER PATH INTEGRAL TECHNIQUE 

In the primitive series approach [c.f. Eq. (8)], the 
"tail" portion of the Brownian bridge is simply discarded. 
Rather than neglecting these terms entirely, it is possible 
to include (approximately) their effects through a num- 
ber of approaches. -One of these is known as the partial 
averaging method. □ Another is a method we term the 
reweighted method introduced in Section IIIB. We note 
that in both methods the n-th order partial sum S" is un- 
changed, its distribution being identical to the primitive 
method one. All methods which preserve the distribu- 
tion of the partial sum S" are referred to by the name of 
the respective series. As such, if the sine-Fourier basis is 



utilized, we shall call the aforementioned approaches the 
PA-FPI and the RW-FPI methods, respectively. 



A. Partial Averaging Method 

Developed initially for the Fourier path integral 
method, the partial averaging technique can be defined 
for all series representations. The key is the indepen- 
dence of the coordinates a^, which physically amounts 
to choosing those representations for which the kinetic 
energy operator is diagonal. Denoting by E„ the aver- 
age over the coefficients beyond the rank n, the partial 
averaging formula reads: 



PPAix,x';/3) 
Pfp{x,x';/3) 




As we mentioned before, the series J^kLn+i'^k-^kiu) is 
again a Gaussian distributed variable of mean zero and 
variance E(i3")^. Using this together with the equal- 



ity it is not difficult to show that formula ( p^ be- 
comes 



Pfp{x,x';f3) 



d^(ai) . . . / d^(a„) exp 



where 



with r^(w) defined by 



exp 



2r2(«) 



V{x + z)dz, 
(16) 



mo 



l-K)-^Afc(7.)2 



fc=i 



(17) 



f3 / Vu,n 

Jo 

I 



V "^0 t^l 



du 



(15) 



There is one property of the partial averaging method, 
of particular note: an application of Jensen's inequalitytHi 
shows that 



p''p\\x,x';P) 
Pfp{x,x';(]) 



d^(ai)... / d^(a„) / d/x(a„+i) exp <j - /3 / Vu,n+ 
dp,{ai) . . 
d/^(ai) . 



xo{u) + J akAk{u) 

V 



d^(a„)exp<^ - (3 d/x(a„+i) / Vu,n+i xoiu) + 

L V "^0 



afeAfc(w) 



d^(a„)exp<{ -/3 / T^„,„ xq{u) + \l ak^kju) 
L V "^0 



fe=i 



du 



du y > 



du 



(18) 



p'^^{x,x';P) 
Pfp{x,x';l3) ■ 
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Therefore, the sequence 

p%^{x,x';(3) <p],a{x.x'-P)<...<pIa{x,x'-(3)<... 

(19) 

is an increasing sequence that converges from below to 
the true density matrix, p{x, x'; /?). 

Let us now consider the problem of choosing the best 
series representation for use within the partial averaging 
framework. We notice that the sequence r^(it), given by 
formula (^, decreases monotonically while pp^(x,x'; (3) 
increases monotonically, as shown by formula (19h. In 
fact, there is a connection between and (^) in the 
sense that the faster the Gaussian spread converges to 
zero, the faster Vu,n{x) converges to the original po- 
tential V{x), and the faster ppj^{x,x';f]) increases to 



p{x,x';P). We note that this observation is general, in- 
dependent of the potential V{x). Of course, one may try 
to optimize pp^(a;, x'; (3) directly, but then the best basis 
will depend upon the potential, an undesirable computa- 
tional feature. We thus conclude that the optimal basis 
for the partial averaging method is the one for which the 
time-average of T^{u) has the fastest decrease to zero, 
i.e. the Wiener or Fourier basis. In this sense, the best 
partial averaging method is the PA-FPI approach. 

We now present one final argument in favor of the 
Wiener basis, an argument that will lead us to a new 
computational approach, the reweightcd FPI technique. 
Remembering the primitive random series method (j^) 
and defining 



Xoc{x,x', a;P) = exp I -(3 V xo{u) + \ akAk{u) 

^ Jo ^ V "^0 



du 



and 



~l3 V xo{u) + \ y^akAk{i 

Jo ^ \ rno 

respectively, we have to the first-order in f3: 

EnXocix,x' ,a; P) - Xn{x, x' ,a; P) « pXn{x,x' ,a; (3) x 




du } , 



.To (u) + W UkAk (u) 

mo 



xq{u) + \ afeAfe(u) 



du. 



(20) 



r 



The rate of convergence for the primitive random series 
method thus depends on the difference between V{x) and 
yu,n{x), which in turn depends on the value of r^j(u). 
Therefore, to a first approximation, the differences be- 
tween the exact and the n-th order FPI density matrices 
depend not on the detailed structure of the respective 
tails, but rather on the spread of the tail series, B"{a), 
a quantity whose time average reaches a minimum for 
the Wiener series. One can readily verify that the term 
of order /3 vanishes for the partial averaging analog of 
formula (^o|), an indication that the technique exactly 
accounts for the extra spread of the paths induced by 
the tail series. 



B. Reweighting Method 

Unlike the partial averaging method, the reweighting 
technique attempts to account for the effects of the tail 
series in a way that does not involve modifying the as- 
sociated potential energy. We shall work out the result 



for the Wiener basis, noting that: (1) the approach can 
be applied to any arbitrary representation, and (2) the 
efficiency of the method will depend upon the specific 
series selected. The basic idea is to replace -B"(a) by an- 
other collection, • • • , 6„), which is supported by an 
n-dimensional underlying probability space. We require 
that: 

1. The variance at the point u of • • • ,6„), de- 
noted by r'„(u), be as close as possible to r^(u). 

2. The variables S'"(ai, • • • , a„) and • • • , 6„) be 
independent and their sum have a joint distribution 
as close to a Brownian bridge as possible. 

One possible candidate for our approach is to choose 
Kibi,--- ,bn) = ^/^j:l=ibA{u), with 6i,---,5„ 
independent identically distributed standard normal 
random variables. Condition 2 above is realized in 
the Ito-Nisio theorem by insuring that the collection 
{cos(fc7ru), a;fe(u)}fe>i is orthogonal, where Wfc(u) is the 



6 



derivative of flk{u). We shall enforce this condition by 
choosing = a„^k sin[(A:+n)7rit] where a„^fc are some 

constants yet to be determined. With the condition 1 
above in mind, and by noticing that in the exact FPI 
representation (Q) the terms of the form sin[(fc + nj)Tru] 
with j > 1 "decouple" as n —^ oo, our intuition tells us 
that a good candidate for is 



-T 



1 



With this choice, the n-th order RW-FPI density matrix 
is given by the formula 



Pfp{x,x';(3) 



dP[a]exp<^ ~ (3 I V 

2n 



xo{u) 



ak(Jn,k sm{kTru] 



k=l 



du 



(21) 



where 



'n,k 



moTT^ 



if 1 < /fc < 



Y.T=o^/ik + inY, iin<k< 271. 



(22) 

The evaluation of the path weights cr^ f. is discussed in 
Appendix D. 

Clearly, our choice of • • • , 6„) is not unique. For 

a better understanding of the quality of the approxima- 
tion, let us compare numerically: 

• r^j(u), the tail variance for the full FPI representa- 
tion and for the PA-FPI method, 

• r'^(u) = fe sin(A:7rM)^, the tail variance 
for the RW-FPI method, and 



{u) — X]fc=n+i '^fc ^i^l^'"'") ' ^^"^ ^^^^ variance 
for the FPI method if it were computed without 
reweighting (i.e. by simply considering the next n 
Fourier terms). 

Fig 1 plots the above variances for n = 9. We notice that 
r^(u) and r'^j(it) are indeed close, much closer than the 
result obtained by simply expanding the primitive FPI 
approach with a similar number of additional terms. 



at each pair of points (x, x'), and their first-order temper- 
ature derivatives converge as fast as ©(l/n"). Generally 
speaking, the aforementioned quantities may have dif- 
ferent asymptotic rates of convergence. However, if the 
potential is smooth enough, our intuition says that this 
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FIG. 1: A plot of the tail variances for the PA-FPI, RW-FPI, 
and 2n-order primitive FPI for n = 9. Notice that the simple 
inclusion of the next 9 terms within the primitive FPI is not 
the optimal strategy. 



is not true. For the case of the harmonic oscillator, we 
shall only verify the convergence of the partition func- 
tion. On the other hand, for numerical simulations it is 
more convenient to compute the average energy of the 
system with the help of the so called T-estimator, which 
can solely be expressed as a functional of the diagonal 
density matrix: 



p{x; P)dx 



(23) 



The above formula can be expressed as the statistical 
average 



/r In <iP[a]Xn{x, a; f3)E'^{x, a; (3) 
/m d^; dP[a]X„(a;, a; /3) 



(24) 



which can be evaluated by Monte Carlo integration. Us- 
ing the notation 



IV. ASYMPTOTIC CONVERGENCE OF THE 
FPI TECHNIQUES 

We say that a given method converges asymptotically 
as 0(l/n") if the partition function, the density matrix 



c„(a, u; /?) = ^ flfco-fe sin(fc7ru) 



fc=i 



to denote the stochastic portions of the current path 
truncated to the first n terms, one easily deduces that 
the T-estimator function for the primitive FPI method is 
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E'^{x,a;l3) = ^ + J V[x + Xnia,u; (3)] du + ^ j V'[x + Xr,{a,u] (3)] Xn{a,u] (3) du, (25) 
while for the PA-FPI method, one obtains 



1 



1 -, rl 



-I 



E^{x,a;(3) = — + / Vu,n[^ + Xn{a,u; P)]du + - / y„ „[a; + x„(a, m; /?)] a;„(a, u; /3) du + 
Jq ^ Jo ' 

1 /■ — 

2 Kni^ + ^nia, u- (3)] Tl{u) du. (26) 

By a simple integration by parts against the coordinate x, one may eliminate the second derivative of the potential 
and obtain the following equivalent PA-FPI energy estimator: 

1 /-^ . . 1 



£'„(x,a;/3) = — + y y„,„[a; + a;„(o, w; /?)] du + - / y„ + a::„(a, m; /?)] a;„(o, u; /3) du + 



't^«,n[2; + 2;„(a,u;/3)]r^(u)du N / V ^,^[x + Xn{a,u; (})\du\ . (27) 



2 

The T-estimator for the RW-FPI method has the same expression as the one for primitive FPI, except for the 
redefinition of the current path 

2n 



c„(a, u; /3) = ^ ak(Tn,k sm^knu). 



k=l 



It is important to note that because of the way we have included the temperature dependence of the path distribution 
in the above analysis, we have obtained direc-tly the so called "virial" forms of the energy estimators. These virial 
expressions have desirable variance propertiest3E3 and are generally preferred for precise Monte Carlo applications. 
The special form ( p7| ) of the PA-FPI energy estimator is numerically advantageous since it does not require the 
evaluation of the second derivatives of the averaged potential. Although we do not study it in this paper because it is 
not a functional of the diagonal density matrix, the H-estimator for the PA-FPI method can similarly be put in the 
simple form: 

E^{x,d;P)^ — + Vix) + -^ / iu-T)^vljx + Xn{a,u;l3)]v',„[x + Xnia,T;P)]dudT. (28) 

The equivalent H-estimator functions for the primitive FPI and RW-FPI approaches look formally the same, except 
that the potential is no longer averaged. The H-estimator is thus properly defined even for potentials that do not 
have second-order derivatives. The reader should notice that the double integral appearing in (^) is really a sum of 
products of monodimensional integrals. We chose this representation for symmetry purposes. The estimator is thus 
the sum of the "classical" energy and a "quantum" correcti on term. 



A. Partition functions for the harmonic oscillator 



In Ref. enough analytical evidence was presented 
to suggest that the asymptotic behavior of the primitive 
and partial averaging FPI methods is controlled at most 
by the values of the second derivatives of the potential. 
Here, we conjecture that this remains true of the RW-FPI 
method, so that an analysis of the harmonic oscillator, 
the simplest potential having a non-vanishing second- 
order derivative, should provide a reliable guess of the 
asymptotic rates for all "smooth" potentials (defined here 
as the potentials having continuous second-order deriva- 
tives). Therefore, we shall study the asymptotic conver- 



gence of the partition function for a one-dimensional par- 
ticle of mass mo — 1 moving in the quadratic potential 
V{x) = x^ 12. We also set ^ = 1 and /? = 1. 

The exact analytical expressions for the harmonic os- 
cillator partition functions are derived in the Appendix B 
for the three methods: primitive FPI, PA-FPI, and RW- 
FPI, respectively. The partition functions of even and 
odd orders have a slightly different convergence behavior 
according to whether 1 — (— 1)" is or 2 (see Appendix B). 
To avoid the appearance of certain oscillations in our 
plots, we shall only compute the odd subsequence for the 
RW-FPI method. Remember, however, that the 2n -I- 1- 
th order RW-FPI approach uses in fact twice as many 
points. To ensure fairness as far as the computational 
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effort is concerned, we shall compare the 2n + 1-th order 
RW-FPI results with those of the 4n + 2-th order primi- 
tive FPI and PA-FPI approaches since, for a given order, 
the former method uses twice as many path variables as 
do the latter techniques. It is convenient to redefine the 
order of the RW-FPI method as being equal to the num- 
ber of random variables used to parameterize the paths, 
in this case: An + 2. In general. 



dP[a]exp<^ -P V 
n I Jo 

2n 



xo{u) 



+ ^ ak<7n,k sin(fc7ru] 



k=l 



du 



(29) 



where an,k is given by formula (^2|) and is evaluated in 
Appendix D. 

Let us assume that we may expand the difference 
Zp^{f3) — Z{f3) as the generalized power series 



Asymptotic rate for the quadratic potential 




FIG. 2: A plot of the indices of convergence for the PA-FPI, 
RW-FPI, and primitive FPI for the quadratic potential. 



Z^M - Z{I3) 



fc=i 



with c ^ 0. For n large enough, it suffices to consider the 
approximation 



Z'-P,{p)-Zip) 



(30) 



By passing to the subsequence An + 2 and taking the 
ratios of consecutive differences, we obtain: 



Zp-\(3)-Z{(3) 



fAn + 2Y l + ci/(4n + 2) 



^Pr 



(/3)-Z(/3) \4n-2j l + ci/(4n-2) 



Next, we take the logarithm and use the approximation 
1/(1 + x) w 1 — x for the last term: 



log 1 



rrAn + 2 



^Pr 



(/3) - Z(/3) 



a log 1 



4n - 2 



log 1 



ci 



4n2 - 1 



We expand the logarithms on the right-hand side of the 
above equation so that the error be of the order 0{l/n^) 
and then multiply the resulting equation by ti^ — 1 /4 to 
obtain 



(n2^1/4)log 1 + 



^Pr 



zp+\p)-z{p) 



na 



-a/2- 



in ■ 



a- 



An-2 

It is convenient to introduce the notation 



ci/4. 



DZ 



Pr 



(/3) = z|,r'(/3)-^r. (/?) 



and set 



(n2 - 1/4) log 1 + 



DZp+\(3) 



n,in+2 
^Pr 



(/3) - Z(/3) 



Since (4n -I- 2)/(4n — 2) « 1 for n large, we conclude 



an - 



a/2- 



21 

4 ' 



(31) 



which shows that ap^ should be asymptotically a straight 
line whose slope gives the convergence order. Here, Z(j5) 
is the exact value of the partition function and the in- 
dex Pr is used to denote the primitive FPI method. Of 
course, similar expressions can be written for the other 
two methods identified by the indices BW and PA. For 
general expressions which apply to any of the techniques, 
we shall use the index Mt. 

Once the asymptotic order is established, we may de- 
termine the value of the constant c by analyzing the slope 
of the equation 



Cm* 



cn - 



c/2 -t- cci. 



(32) 



where 



cl,, = (4n -f 2.r(n + 1/2) [Zt?+^(/3) - Z{(i)\ 



The asymptotic behavior implied by (^2|) can easily be 
established by replacing n by 4n -I- 2 in equation ( |30| ) . 

Fig. 2 shows that the linear region predicted by our 
analysis is quite rapidly reached for the harmonic os- 
cillator. One easily notice that the PA-FPI and RW- 
FPI methods have similar asymptotic behavior, while the 
primitive FPI approach has a slower rate of convergence. 

The asymptotic slopes are computed as the slope of 
the line that best fits the last [A^/3] values, where N is 
the number of data points calculated. We assume that 
we computed enough points so that the last \N j 3] are in 
the asymptotic region. Euler least-square fit gives then 
the value 



OLMt 



(33) 
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where the summation is done over the last data 
points. Of course, the exact value for a is the limit as 
— > oo of the right-hand side of the above formula. 
For N = 12, (H) gives: apr = 1.002, apA = 3.007 and 
o^RW = 3.008, suggesting that the asymptotic behavior is 
0{l/n) for the first, and 0{l/n^) for the last two meth- 
ods, respectively. 

The constants c are calculated in a similar fashion and 
the numerical values for iV = 12 are: cpr = 0.049, 
cpA = —7.933 • 10"'^, and cpw = 0.887, respectively. 
Therefore, the partial averaging method is superior to 
the reweighted method in the sense that it has a smaller 
convergence constant (smaller in modulus). 

Theoretically, if we can compute the difference between 
successive values of the partition function with sufficient 
precision, we can improve the convergence of any of the 
FPI methods by using better estimators. For first-order, 
the result can be obtained as follows: formula ( |30| ) shows 
that 



Z 



4n+2 
Mt 



(4n + 2)" 



(34) 



converges to the exact answer as fast as 0(1/ (4n+2)"+^) 
and therefore, the last equation is a better estimator as 
far as the asymptotic behavior is concerned. Given that 
the convergence exponent a is known, the constant c can 
be approximately (but arbitrarily exactly as n — > oo) 
evaluated from the equation: 



n nrAn+2 



(P) 



(4n-h2)"- (4n-2)" _ c a 
(4n-h2)"(4n- 2)" (4n + 2)" n 



Solving for c and replacing in ( [3^ ) , one ends up with the 
first-order corrected estimator 



P74ti-(-2 



Zil+^{(3) - ^DZ'^^+\(3) (35) 



The second-order estimator can be derived by apply- 
ing the first-order correction to the first-order estimator. 
One easily computes: 



(2a + l)n 4n+2(a-. 
a{a + l) 



-DZt 



a{a + 1) 



a(a + l)""^^* 
'DZt}r\P)-DZtit\P) (36) 



The asymptotic convergence of this estimator is 
C'(l/n"+2). 

In principle, one can continue this process beyond 
second-order. However, as we shall see in the Ap- 
pendix E, such higher order estimators are of little practi- 
cal value. To demonstrate the behavior of the corrected 
estimators, we compute the convergence exponents for 
the primitive FPI using the corresponding analog of equa- 
tion (^l|). Fig. 3 clearly shows the difference in the rate 
of convergence for the original and corrected estimators. 



Asymptotic rates for different estimators 
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FIG. 3: A plot of the exponents of convergence for the three 
^-estimators. The method employed is primitive FPI as ap- 
plied to the quadratic potential. 



The numerical values are az = 1.002, apz — 1.997, and 
ctgz = 2.958, demonstrating our predictions. From now 
on, we shall refer to the original, unaccelerated estimator 
as the zero-order estimator. 



B. A numerical example: The quartic potential 

As we said in the beginning of this section, for numer- 
ical purposes it is convenient to study the convergence of 
the T-method energy estimator in the virial form, which 
can be computed by Monte Carlo integration. As we 
explain below, the numerical study of the asymptotic be- 
havior is not a computationally easy task, especially for 
those methods that have rapid asymptotic convergence. 
More explicitly, let us take a look at the following analog 
of (H): 



'■Mt 



OLMtn - OLMtl'^ ■ 



Cl,Mt 



(37) 



where 



Mt 



(n^ - 1/4) log 1 + 



E' 



4n-2 
Mt 



E 



4n+2 
Mt 



^Mt - ^ 



For the partial averaging method, we suggested that the 
difference E'^^'^ — E decays to zero as fast as l/m?. In 
turn, the differences E^^'^ — Ep'^'^ between consecu- 
tive terms decay to zero as fast as 1/n^. It is thus clear 
that faster rates of convergence of the method require 
greater precision in the evaluation of the terms E'^^^ . 
If we assume an independent sampling of the probability 
density shown in formula (p^), the error in the Monte 
Carlo evaluation of Ep^^ is AEp/^/VN, where N is 
the number of Monte Carlo sampling points and Ai?p^^ 
is the standard deviation. This error should satisfy the 
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inequality: 



\AEp/^\/VN«\Ep^- 



77i4n+2| 



Asymptotic rate for the quartic potential 



It follows that the number of Monte Carlo points nec- 
essary to insure a given relative error for apA scales at 
least as badly as iV oc as a function of the number 
of Fourier coefficients. The same is true for RW-FPI, 
while for the primitive FPI we only need N (x n'^. We 
emphasize that this scaling is related to our immediate 
task of establishing the asymptotic rates of convergence 
and is not an issue that would arise in typical numerical 
applications. 

The second observation we make is that the ratio 



\E 



Mt 



Et-+'\/\AEZ\+'\ 



increases as the temperature is dropped. Consequently, 
we would like to conduct our model computations at low 
temperature, where the quantum effects are big enough 
so that the differences between consecutive terms are sig- 
nificant. At high temperature, the classical limit is a 
good approximation and these differences may be smaller 
than the statistical errors we are able to achieve. We are 
therefore forced to conduct our computations in the "un- 
favorable" range of temperatures, and in general, we need 
to study groundstate problems. 

We hope this is enough rationale to justify the need for 
a special Monte Carlo integration scheme capable of accu- 
rately sampling the low temperature distributions with 
good efficiency and low correlation, at least for certain 
classes of simpler systems. One such scheme is discussed 
in Appendix E, and it generally applies to the class of 
single-well potentials. 

For comparison purposes, we shall also compute the 
T-estimator energies for the trapezoidal Trotter method. 
Expressions similar to those presentejd here for the FPI 
methods were deduced by Coalsontil and employed by 
Mielke and Truhlaifl as the TT-FPI method. We shall 
keep this name in the present paper, though, as defined 
here, the TT-FPI approach is not an FPI method because 
the n-th order partial sum 5*" is not the one for the prim- 
itive FPI. The importance of this method consists of the 
fact that its asymptotic rate of convergence is 0{l/n'^) 
for smooth enough potentials, being the fastest primitijie 
method to date that leaves the potential unchanged.Ej 
We do not present this scheme in the present paper and 
for further information we refer the reader to the cited 
literature. 

The prototype system studied in this work is the quar- 
tic potential V{x) = We set h = 1 and toq = 1 and 
(3 — 10. The groundstate of the quartic potential was 
evaluated by variational methods to be Eq = 0.530181, 
while the average energy at the temperature correspond- 
ing to /3 = 10 is i? = 0.530183. We computed the average 
energy for the sequence An + 2 with 1 < n < 12, corre- 
sponding to the actual numbers of Fourier coefficients 
6, 10, . . . , 50. In these calculations, the number of points 
employed in the Gauss-Legendre quadrature scheme was 




FIG. 4: The straight lines drawn represent the linear least 
square fit for the last four data. Their slopes give the conver- 
gence exponents for each method. 



200. We used 1.25 • 10^ Monte Carlo points for primitive 
FPI, 2.5 • 10^ for TT-FPI, 5 • 10^ points for RW-FPI, and 
2 • 10^ points for PA-FPI calculations, respectively. The 
values of i?p^^ were previously computed in a quarter of 
these numbers during a "warm-up" period, but we con- 
tinued to improve them during the main Monte Carlo 
procedure. Table || of Appendix F summarizes the re- 
sults of the computer evaluations. The differences be- 
tween successive energ y te rms were computed with the 
help of the estimator (E9). T he errors were computed 
with the help of the formulae ( E3 ) for the average ener- 
gies, and (Ell) for the estimated differences. 

Fig 4 shows the behavior of the functions a^,j^ for the 
four methods. Among the non-averaged methods, we re- 
mark that the primitive FPI approach reaches its asymp- 
totic behavior faster than the TT-FPI method, which in 
turn reaches its asymptotic region faster than the RW- 
FPI technique. This behavior is shown in Fig 5, which 



plots the current slope d^j^ — a^,^ . Although the RW- 
FPI method did not reach its final asymptotic behav- 
ior, the trend is clear. The computed convergence expo- 
nents using the last four data points are: apA = 3.082, 
apw = 2.917, aTT = 2.071, and apr = 1.019. There- 
fore, we conclude that the asymptotic convergence of the 
methods is 0{l/n^) for PA-FPI and RW-FPI, ©(l/n^) 
for TT-FPI, and 0{l/n) for primitive FPI. Lastly, it 
is worth comparing the convergence constants for the 
PA-FPI and RW-FPI methods since they have the same 
asymptotic convergence exponent. The numerical values 
are cpA = 59.4 and cpw = —736.7, showing that the 
PA-FPI method is over 10 times faster than the RW- 
FPI method. This is in agreement with the observations 
made for the partition function of the harmonic oscillator 
in the previous section. The PA speed-up of the conver- 
gence is important, especially with respect to minimizing 
the number of path variables required in practical appli- 
cations. 
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FIG. 5: The current slopes for each method should ideally 
converge to 3 for PA-FPI and RW-FPl, 2 for TT-FPI, and 1 
for primitive FPL 



it is the best way of improving the asymptotic behavior of 
the primitive FPl method (at the cost of computing the 
Gaussian transform of the potential). The TT-FPI and 
RW-FPl methods increase the order of convergence of 
the primitive FPl to ©(l/rt^), and 0{l/n^), respectively, 
without increasing the variance of the corresponding esti- 
mators. It should be noted that, unlike the complete par- 
tial averaging approach, the reweighting method does not 
require the Gaussian transform of the potential. While a 
final decision awaits detailed future studies, we anticipate 
that this latter feature of the reweighting approach will 
be beneficial for applications where the Gaussian trans- 
form is either formally ill-posed and/or computationally 
difficult to obtain. Finally, as discussed in Appendix E, 
the first and the second-order estimators also improve 
the asymptotic convergence. Although both have larger 
variances, the first order estimator appears computation- 
ally feasible since its variance decreases with the number 
of Fourier coefficients. 



CONCLUSIONS 



In this paper, we have shown that the best series rep- 
resentation (with respect to asymptotic convergence) for 
use in Monte Carlo path integral methods is the Wiener 
sine-Fourier series. Both the RW-FPI and TT-FPI meth- 
ods are not series representations and we suggest that the 
latter also falls in the category of reweighting techniques. 
The partial averaging technique has the asymptotic con- 
vergence 0{l/n^), with a small convergence constant and 
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APPENDIX A: ITO-NISIO THEOREM 

Theorem 1 (Ito-Nisioi) Let {Afc(T)}fc>o he any or- 
thonormal basis in L^[0, 1] such that Ao(r) 1, let 



is uniformly convergent almost surely and equal in distri- 
bution with a standard Brownian bridge. 



AfeH 



Afe(T)dr, 



and let a = {ak}k>i be a sequence of independent iden- 
tically distributed (i.i.d.) standard normal random vari- 
ables. Then, the series 



APPENDIX B: HARMONIC OSCILLATOR 



J 



The 2rt-th order primitive FPI approximation of the partition function for an harmonic oscillator centered at the 
origin has the expression 



Zjr^^ f dx f dai... / da2nP"pr{x,a;(3), 



where 



afefTfc sin(fc7rM) 



k=l 



du 



By explicitly computing the integral over t and then completing the square, one obtains 



p|,"(x,a;/3) 



(27r)-"exp 



2n 



k=l 



£,kak 



kn 



exp ■ 



1-E 

k=l 



Pmouj'^al /I - (-l)*^ 



i 



kn 



(Bl) 



where 

ek = l + l3moiu'al/2. 

For use as a trial density in Monte Carlo simulations, it is convenient to replace the last factor by its limit n — > oo: 



exp 



— ^ — tanh 



huj 



1 



(27r)' 



■ exp • 



2n 

-Y 

2 ^ 

k=l 



ik^k 



f3moUj'^x 1 - (-l)* 



kn 



(B2) 



It is not difficult to show that the trial densities for the primitive FPI and PA-FPI methods are identical (after 
normalization) but we shall employ formula (B2) for the R W-F PI technique too. Practice shows that the penalty for 
considering the last two approximations is minimal, while (B2) has some advantages with regard to the organization 
of the computations. 



To evaluate the partition hmctions for the primitive FPI approach, we integrate (Bl) and obtain 

, -1/2 



1 



2n 



.k=l 



2n 



1-E 



f3mouj'ai fl- (-1) 



fe=i 



kn 



We leave for the reader the simple task of showing that the 2n-th order PA-FPI density matrix has the form 

2 2n 



/ n' 



27r2 



6 



k=l 



fc2 



(B3) 



(B4) 
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The RW-FPI method's partition function is similar to the one for the primitive FPI method and is given by 



nrr 



2n 



1-E 



— 1 ^n,i 



kn 



-1/2 



(B5) 



where 



I 



APPENDIX C: METROPOLIS SAMPLING 

The Metropohs et al. samphng of a general prob- 
ability density p{x) with x €: fl consists of generating a 
homogeneous Markov chain having the transition proba- 
bility density 

t{x'\x) = A{x'\x)T{x'\x) + 
6{x' -x) [ [l-A{y\x)]T{y\x)dy, (CI) 



where T(x'\x) is a trial transition probability density 
which would generate an irreducible chain by itself, 
d{x' — x) is the Dirac function, and the acceptance prob- 
ability A{x'\x) is given by the formula 



= min < 1, 



Pi^')T{x\x') 
p{x)T{x'\x) 



This choice of ^(a;'|a;) is one of the many possible which 
satisfy the condition 

A{x'\x)T{x'\x)p{x) = T{x\x')A{x\x')p{x'). 

The last relation implies that the Markov chain of transi- 
tion probability density t{x\x') satisfies the detailed bal- 
ance condition 

t{x'\x)p{x) = t{x\x')p{x'), 

which by integration against x' and use of the normal- 
ization condition J^^ p{x)dx = 1 shows that p{x) is a 
stationary distribution of the transition kernel t{x'\x). 
Moreover, it can be shown that the associated Markov 
chain is ergodic and that this—implies that p(x) is the 
unique stationary distribution.tfl Let us consider the sta- 
tionary sequence Xq ,Xi,... with Xq having the distribu- 
tion density p{x) and Xn having the conditional density 
P{Xn = x'\Xn-i — x) = t{x\x'). One can generate 
a sample xo,xi^ . . . starting with any point xq, by the 
Metropolis algorithm: 

1. given Xn, generate Xn+i from the probability den- 
sity T{x\Xn); 

2. compute A{xn+i\xn); 

3. generate a random number q uniformly on [0, 1]; 



4. if q < A{xn+i\xn)i accept the move; otherwise, re- 
ject it. 

For the expected value ]E(/) = p{x)f{x)dx, Birkhoff's 
ergodic theorem (Theorem 2.1, Chapter 6 of Ref. ^) guar- 
antees that 



n— 1 

-5:/(x.)-E(/) 



(C2) 



fe=0 



almost surely. In words, the probability that we may gen- 
erate a sequence xo,xi, . . . by the Metropolis algorithm 
such that 



n — 1 

J2fix,)^E{f) 



k=0 

is zero. In fact, if the variance of f{x) is finite 

ao^(/)=E(/-E/)2<oo, 

a central limit theorem holds. Since the random vari- 
ables f{Xo) and f{Xn) have the same distribution, their 
correlation coefficient takes the form 



rnif) 



E[f{Xo)f{Xn)]~E{ff 



Explicitly, let us introduce the notation 

..V|.)=/d.,.../d..._,.,.'|.,)....(..-W. 

Jn Jn 

with T°{x'\x) = 6{x' - x) and t^{x'\x) = t{x'\x). Then, 



E[/(Xo)/(X„)] = / dx / dx'pix)T"ix'\x)f{x)f{x'). 
Jn Jn 

In practice, we can evaluate these expectations, and 
therefore the correlation coefficients, again with the help 
of Birkhoff's theorem: 

^ fe-i 

E [/(Xo)/(X„)] = hm - J{x,)f{x,+n). (C3) 



J=0 



In these conditions, it can be shown (Theorem 7.6, 
Chapter 7 of Ref. |) that 



(C4) 
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where ^ has the standard normal distribution and 



^'(/)=<^o(/) 



.(/) 



(C5) 



If the samphng were independent, the correlation coef- 
ficients would vanish and we would recover the classical 
central limit theorem. In practice however, the correla- 
tion coefficients are positive, many times having a slow 
decay to zero and the independent sampling may be con- 
sidered a fortunate case. Without entering the details, 
we mention that there are two factors that contribute to 
large correlation coefficients: a) a strongly correlated pro- 



posal T{x'\x) and b) a low overall efficiency. The overall 
efficiency (or the acceptance ratio) is defined as 



dx' dxpix)Aix'\x)T{x'\x) (C6) 
SI Jq 



and represents the fraction of moves accepted. Therefore, 
if the overall efficiency has large enough values {Ac > 
0.2), it is a good idea to use an independent proposal 
from a trial probability pt^{x). If ptri^) ~ p{x) and 
J{x) is smooth enough, we may approximately relate the 
correlation coefficients to the overall efficiency as follows: 
from the relation (CI), we easily compute 



TlU) - 1 



j^dx'j,,dx[f{xf - f{x)f(x')]p{x')ptr{x)A{x'\x) 

dx% dxifixr - f{x)f{x')]p{x')p{x) 



r 



where 



A(.'|.).min|l,4:^^ 

I P{x)ptr[x') 



(C7) 



Using the approximation p{x')ptr{x)A{x'\x) « 
Ac p{x)p{x' )^ the right-hand side simplifies to 
ri (/) w 1 — Ac. In general, by a similar line of 
thought, one may argue that r„(/) « (1 — ^c)". The 
formula (C5) takes the approximate value 



^'(/)«^o(/) 



1 + 2Y^{1 - Acf 



k=l 



(C8) 

Therefore, the bigger the acceptance probability, the 
faster the convergence of the Monte Carlo procedure. In 
the limit Ac = 1, we recover the independent sampling, 
but a quick look at formula (C7) shows that in this case 
Ptr{x) = p{x). 



APPENDIX D: COMPUTATION OF THE PATH 
WEIGHTS gI FOR THE RW-FPI METHOD. 



Hurwitz ^-function, usually implemented by many math- 
ematical libraries. Alternatively, h(a;) can be evaluated 
via the trivial identity 



\i{x) = C(2)-2a;C(3)-|-3a;2C(4)-a;^ 



3 



4j 3a; 

\ (JTWf 



(D2) 



where C(s) is the Riemann ^-function 



C(5) 



oo 



We have C(2) = C(3) ~ 1.2020569031596, and 

C(4) = 7r'*/90, with the last series in ( p2| ) converging 
quite fast. More precisely, the error in the evaluation 
of h(x) committed by truncating the series to the first n 
terms is easily seen to be smaller than J2j>n ^/j^ — 1/"'' 
uniformly on the whole interval [0, 1], so that summation 
over the first 100 terms gives the value of h(x) with an 
error of at most 10~^. This error is sufficiently small for 
our applications. 



If n < fc < 2n, we have 



2(3K 



'n,k 



where 



2 oo 



2/3^2 I fk~n 
h 



(Dl) 



^(-) = g(7T^- 

Clearly, the values of the function h(a;) are only needed 
over the interval [0, 1] and they can be evaluated via the 



APPENDIX E: A SPECIALIZED MONTE CARLO 
SCHEME 

As suggested in Appendix C, the use of an indepen- 
dent trial distribution in the Metropolis algorithm is a 
good strategy provided that we are able to find a good 
approximation pj"+2(a;, a; /3) to the density we need to 
sample in this case, 



p^^+\x,a;P) 



(2^) 



x:r'(x,a,/?) 



2n+l 



■ exp 



4n+2 
k=l 



(El) 
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The quartic potential 















rriQU^ 2 
2 ^ 










X 



FIG. 6: A plot of the quartic potential (solid line) and its 
best variational quadratic approximation. Here, mo = 1 and 
= 1.442. 



This approximation may be taken to be the similar ex- 
pression for a harmonic oscillator potential mQijj'^{x — 
A)"^ /2, because we know how to generate an independent 
sample of this. In order for the approximation to work 
well for many of the single well potentials of interest, we 
optimize the parameters w and A to obtain a best fit in 
the sense of increasing the overall acceptance ratio. How- 
ever, since we are analyzing groundstate problems, suf- 
ficiently good approximations can be obtained from the 
Ritz variational principle. Thus, we look for the parame- 
ters Lo and A which realize the minimum of the functional 



where 



fmoU}\ 



1/4 



exp 



2h 



{x-AY^ 



is the groundstate eigenfunction of the trial harmonic 
potential and 



H = 



2to( 



■A + V{x) 



is the Hamiltonian of the original single well potential. 
By a translation of the reference system, we may assume 
that the optimizing parameter A is zero. For the case of 
the quartic potential V{x) = x'^/2, the best optimizing 
parameters are uj = 1.442 and A = 0. Fig. 6 plots the 
quartic potential and its best quadratic approximation. 

Rather than using the 4n -I- 2-th order probability den- 
sity of the best harmonic reference as the trial density, 
it is more convenient to use the slightly modified for- 
mula (B2) of Appendix B. The advantage is that (B2) 
is the exponential of a series. As such, if we generate 
the vector {x,ao, ...,a4„+2) from the probability density 
a; P) given by (p2|), we can use the vectors of the 



form (x, ao, 04^+2) with k < n for the paths of smaller 
length because it is clear that these vectors are drawn 
from the distribution pl^^'^{x, a; (3). The time saved with 
the generation of random numbers fully compens ates the 
slight decrease in the acceptance ratio. We use (B2) for 
all FPI methods in the following examples. 

For PA-FPI and primitive FPI, there is another advan- 
tage in using the trial density (B2). A large portion of 
the computational time is spent with the construction of 
the paths 



in+2 



i+2(w;/3) = ^ ak(7kSm{kTTu), 



k=l 



especially for large n. However, if the trial probability 
density (B2) is used, we can employ the recurrence for- 
mula 



4TI+2 

a4„+2('«;/3) = a4ri-2 (-";/?) -I- ^ 0^0-^ sin(fc7ru). 

fc=4n-l 

Therefore, the time necessary to construct all the paths of 
length 4:k + 2 with 1 < fc < n at a given point t scales like 
0{n) instead of 0{n^). This is especially important for 
the PA-FPI method, which has the fastest convergence 
and for which a large number of Monte Carlo steps is 
necessary to establish the asymptotic convergence rate. 
Unfortunately, since the paths for RW-FPI are not series, 
we cannot employ the same strategy for this method. 

As shown in the Appendix C, the advantage of our 
Monte Carlo strategy consists of the fact that it has low 
correlation provided that the acceptance ratio is large. 
To a first approximation, the statistical error in the esti- 
mation of the energy is [we employ the usual 2a defini- 
tion for the error, corresponding to a confidence interval 
of 95.4%1 



Err. 



2(Jo(E 



An+2\ 
Mt ) 



N 



2 

'Ac 



1 



1/2 



(E2) 



where N is the number of Monte Carlo points, ctq {E'^^'^^ 
is the variance of the T -est imator function, and Ac is the 
acceptance ratio [see (|C8| )]. A more precise formula is 
given by (F" 



Err. {Etir) 



2ao{E 



Mt 



N 



-1 1/2 



2Y.r. 

k=l 



(E3) 

and we have shown in Appendix C how the correlation 
coefficients can be evaluated during the Monte Carlo pro- 
cedure. However, we can use (E2) to find the number of 
steps after which the correlation becomes negligible. In 
the case of the quartic potential, the acceptance ratio 
was bigger than 0.6 for all simulations performed. Since 
2/0.6 - 1 = 2.333 « 1 + 2^^i 0.4*^ = 2.332, we may 
safely truncate the series in (|E3D to the first eight corre- 
lation coefficients and we shall do so for all computations 
concerning the quartic potential. 
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Another important aspect in our computations is the 
numerical evaluation of the one-dimensional time aver- 
ages that are iniialved. This issue was extensively studied 
by Sabo et. al.j^ who concluded that a Gauss-Legendre 
quadrature in a number of points equal to three times the 
number of Fourier coefficients should suffice for most ap- 
plications. We also employ the Gauss-Legendre quadra- 
ture scheme, but in a number of points equal to four 
times the maximum number of Fourier coefficients com- 
puted. Extensive computer observations show that the 
relative error in the evaluation of the T-estimator func- 



tion is smaller than 10~® for the quartic potential. Of 
course, for real-life applications we do not need such a 
precision but here it is important to rule out any factor 
likely to alter the asymptotic law of convergence. 

Earlier in this section, we saw that the scaling of the 
number of Monte Carlo points with the number of Fourier 
coefficients was dictated by the decay of the differences 
^Mt"^ - £^Mt"^^ which we shall denote by DE*^+^. We 
shall improve on this fact by directly evaluating these 
differences with the help of a biased estimator. Define 



and 



Next, define 



r>4n+2 



ri^+'{x, a; (3) = X^i^^x, a; f3)/Xl^+'{x, a; (3) 
Im In ^P["]^'in+2{x, a; /9)r^"+^(a;, a; (3) 



Ir^^ /odP[a]X4„+2(a;,a;/3) 



DElf^^ix,-a;f3) = i?i„!"2'(a;,a;/3)rl:?+^(x,a;/3)/i?lJ+^ - Ei;,X^{x,a; P). 



r,T,Mtr 



j4n+2 j-,T,Mt/ 



It is a simple exercise to show that 



E 



Mt 



Im In dP[a]X 4n+2{x, a; (3)DE^;^j^2 (^^ P) 
Jjjdx ^^AP[a]Xin+2{x,a\ (3) 



(E4) 



(E5) 



(E6) 



(E7) 



A biased estimator for the function ( p^ ) can be constructed as follows: assume you are given a sequence {xk,ak) 
with 1 < fc < N , which samples the probability distribution (El). At step k, compute 



k 

^m"^' = iE^1«'(^.'%-;/3) and Err,(fc,i?t?+2). 



and construct the function 



Then, the biased estimator is defined by the well-known recurrence formula 

\k - l)D£;t7/''"+' + DE:i:^^{x,,au^P)\ /h 



n j7.fc,4n+2 



(E8) 
(E9) 



starting with DE^f^"''^'^ = 0. Clearly, Di?^^*"^^ converges to DE"^^"^ as k gets large. 

The bias in (E8) is due to the fact that we do not use the exact value of R'^^'^ but its unbiased statistical estimator. 
However, for large enough k, it is not difficult to justify the estimate: 



,T,Mt, 



It follows then that the error due to bias is at most 



\El;^'f^{x, a; ri^+\x, a; (3) Err,(fc, 



< 



R 



An+2 
Mt 



R 



4n+2 
Mt 



Eri;{N,DEtl+') 



]_^\j^^^Xk,ak;f3)\rl^+\xk,ak;P) Err,(fc, 



k=l 



R 



k,4n+2 
Mt 



R 



,fe,4n+2 
Mt 



(ElO) 



The total error is then obtained by also adding the statistical error computed with the help of the formula (E3): 



Err(iV, DEl';+^) = Err, (TV, DE 



Mt 



Ern{N, DE 



4n+2\ 
Mt I 



(Ell) 
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In the present paper, we pre-computed a start value of R^^"^ using a quarter of the number of Monte Carlo points 
during the warm-up step and then continued to improve the value in the main procedure. In these conditions, one 



may argue that the error for the difference (E7) satisfies the inequality 



where 



E(|i?: 



4n-2| 
Mt I 



/r /n dP[Q]X4n_2(a:, a; /?) |g4„_2 {x, a; (3) 
/r In <iP[a]X4n-2{x, a; (3) 



Formula (E12) helps us explain why the use of the 
biased estimator ( |E^ ) is advantageous. Had we directly 
evaluated the difference 



DE 



4n+2 
Mt 



= E' 



4n-2 
Mt 



E 



the error would have been 



Err{N,DEt-+') 



Err, 



An+2 
Mt ' 



Err, 



(E13) 



{N,E^^+^ 



Mt ) ■ 

(E14) 

Notice however that both r^"j^^(a;, a; /?) and 

DEj^^2ix,<i', P) converge to 1 and 0, respectively 
as n ^ oo. In turn, their variances (which control 
the statistical errors) converge to zero. Clearly, this is 
not the case for the variance of the T-method energy 
estimator. More precisely. Table || presents strong 
numerical evidence suggesting that the decay of their 
standard deviations is as fast as 0{l/n^) and we expect 
this to be true for all smooth enough potentials. This 
implies that for a fixed but large number of Monte Carlo 
points N, the error in (Ell) has the asymptotic behavior 



This readily implies that the use of the estimators of 
order one and two does not change the scaling of the 
number of Monte Carlo points needed to achieve a given 
error threshold for the estimated energy with the number 
of Fourier coefficients. The net result is an improvement 
in the asymptotic behavior for the estimators of order 
one and two. However, in the case of the second-order 
estimator, we notice an increase in the variance of the es- 
timator which may be quite large for practical purposes. 
For the first-order estimator there is no asymptotic in- 
crease in the variance, which makes it more suitable for 
practical applications. In fact, the first-order estimator 
may also be used for potentials that do not have con- 
tinuous second-order derivatives but for which the decay 
with the number of Fourier coefficients implied by ( E15| ) 
can be replaced by the slower one 



V fAT mpT.Mt. _ const 



N 



Err(iV,i^ii;r4^*; 



const 



N 



(E15) 



The importance of (E15) is twofold. First, it shows 
that if the estimator (E9) is used, the scaling of the num- 
ber of Monte Carlo samples with respect to the number 
of Fourier coefficients is now determined by the decay of 
-^Mt^^ — i? to zero. More precisely, we have N (x for 
PA-FPI and RW-FPI, N n'^ for TT-FPI, and N 
for primitive FPL 

Second, the errors of the estimators of order one and 
two [see (^H]) and (^)] have the asymptotic behavior: 



Err(^,i^i?J„%*) 



= Err 



a n\/N 



Err(7V,i?X;%*) 



(E16) 



and, 



Err(7V,5£;J„*^2*) 
1 



const 



a{a + 1)VN 



2a + I 



Err,(iV,i?J„^^2*) + 



Err, (TV, i;, 



n-1 

T.Mt-, 



2 + 



^4n+2 

2^2 



(n- 1)2 
4 • const 



4n+2 ) 



cx{oi + l)ViV 



(E17) 



Finally, in the cases where it cannot be utilized as an en- 
ergy estimator because of an unduly large variance, the 
correction term brought in by the first-order estimator is 
still useful as a measure of how far the zero-order esti- 
mator is from the true result. 

The reader may work out the expression for the esti- 
mator of order three and see that in this case the scaling 
is changed. This explains our earlier assertion that the 
estimators of order three or more are of little practical 
value. 



APPENDIX F: 



TABLES OF NUMERICAL 
VALUES 



The following tables contain the numerical results de- 
scribed in Section IVB. See that discussion for the details. 
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TABLE I: Average energies, estimated differences, and their statistical error for the quartic potential at /3 = 10. The variational 
energy is 0.530183. 



n 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


12 


Average energies 




0.302878 


0.365234 


0.401528 


0.425003 


0.441342 


0.453379 


0.462581 


0.469834 


0.475704 


0.480548 


0.484613 


0.488071 




0.343731 


0.416263 


0.454541 


0.476808 


0.490728 


0.499978 


0.506376 


0.510972 


0.514391 


0.516994 


0.518994 


0.520602 


p4n+2 
^RW 


0.351676 


0.432846 


0.473011 


0.493918 


0.505667 


0.512786 


0.517363 


0.520451 


0.522627 


0.524201 


0.525370 


0.526247 


p4n+2 


0.596947 


0.552843 


0.541042 


0.536268 


0.533916 


0.532629 


0.531862 


0.531383 


0.531069 


0.530854 


0.530701 


0.530593 


Enti.nutt<:d (Itjj'crcii.ccs 




-.124981 


-.062354 


-.036316 


-.023475 


-.016353 


-.012025 


-.009205 


-.007265 


-.005872 


-.004848 


-.004062 


-.003452 




-.147346 


-.072843 


-.038307 


-.022241 


-.013929 


-.009222 


-.006401 


-.0046052 


-.003424 


-.002593 


-.002013 


-.001589 




-.162238 


-.080976 


-.040075 


-.020896 


-.011746 


-.007111 


-.004573 


-.003100 


-.002176 


-.001575 


-.001168 


-.000886 




0.549021 


0.044095 


0.011781 


0.004772 


0.002347 


0.001285 


0.000763 


0.000481 


0.000316 


0.000216 


0.000151 


0.000109 


Statistical errors for energies{2a) 


^Pr 


0.000088 


0.000084 


0.000081 


0.000080 


0.000078 


0.000078 


0.000077 


0.000077 


0.000076 


0.000076 


0.000076 


0.000076 


j^4n+2 


0.000056 


0.000057 


0.000056 


0.000055 


0.000054 


0.000054 


0.000054 


0.000053 


0.000053 


0.000053 


0.000053 


0.000053 


p4n+2 
^RW 


0.000043 


0.000042 


0.000041 


0.000040 


0.000039 


0.000038 


0.000038 


0.000038 


0.000038 


0.000037 


0.000037 


0.000037 




0.000024 


0.000023 


0.000021 


0.000020 


0.000020 


0.000020 


0.000020 


0.000019 


0.000019 


0.000019 


0.000018 


0.000018 


Statistical errors for differences {2a) 




0.032356 


0.000304 


0.000095 


0.000055 


0.000036 


0.000026 


0.000020 


0.000016 


0.000012 


0.000010 


0.000009 


0.000007 


7-> j-An+2 

UJZjrprp 


0.054577 


0.001325 


0.000199 


0.000092 


0.000061 


0.000046 


0.000036 


0.000030 


0.000025 


0.000021 


0.000018 


0.000015 


n Ei4n + 2 
^^RW 


0.037799 


0.001600 


0.000333 


0.000074 


0.000045 


0.000032 


0.000025 


0.000019 


0.000016 


0.000013 


0.000011 


0.000009 




0.004595 


0.000080 


0.000027 


0.000014 


0.000009 


0.000007 


0.000005 


0.000004 


0.000003 


0.000003 


0.000002 


0.000002 



TABLE II: Standard deviations for r'^^^{x, a; j3) and DE'^^1{x, a; 13) and their asymptotic convergence exponents a. 



n 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


12 


a 


Primitive FPI 


4n+2 
Pr 


1.771976 


0.437892 


0.226571 


0.143267 


0.099735 


0.073744 


0.056769 


0.045052 


0.036617 


0.030341 


0.025542 


2.027 




0.907676 


0.222962 


0.116896 


0.075168 


0.052794 


0.039381 


0.030502 


0.024296 


0.019827 


0.016464 


0.013907 


1.985 


TT-FPI 


4n + 2 
' TT 


9.411364 


1.149898 


0.499554 


0.326273 


0.241282 


0.187604 


0.151356 


0.125104 


0.105345 


0.090015 


0.077867 


1.822 




5.259389 


0.856299 


0.342638 


0.204763 


0.148354 


0.112764 


0.089384 


0.072996 


0.060898 


0.051628 


0.044358 


1.825 


RW-FPI 


„4n + 2 
' RW 


17.63636 


1.953891 


0.543819 


0.324225 


0.228415 


0.171841 


0.134647 


0.108553 


0.089442 


0.074985 


0.063776 


1.961 




8.421340 


2.378073 


0.393884 


0.210005 


0.142382 


0.105825 


0.082139 


0.065720 


0.053826 


0.044884 


0.037983 


1.999 


PA-FPl 


„4n-|-2 
' PA 


1.164360 


0.309710 


0.182490 


0.122904 


0.088704 


0.067010 


0.052355 


0.041982 


0.034380 


0.028648 


0.024221 


2.100 


j^pT,PA 


1.596240 


0.291506 


0.142120 


0.087909 


0.060274 


0.044039 


0.033599 


0.026465 


0.021387 


0.017629 


0.014785 


2.013 



